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超级 计算 面临 的 挑战 


及 其 对 未 来 数值 算法 设计 的 可 


摘要 超级 计算 机 浮 点 运算 能 力 的 迅速 提高 


USC 


HR 


沟 ”。“ 浮 点 效率 鸿沟 ”反映 了 传统 数值 算法 与 


处 理 器 上 ， 
这 一 轮 技 术 变 革 带 


中 


来 的 新 的 免费 午餐 " 


角度 出 发 ， 发展 一 种 新 的 高 效 且 高 精度 无 额 
法 (以 及 有 限 差分 ) 等 一 类 新 兴 数 值 方法 中 所 


有 


精度 “双赢 ”的 新 型 计算 模式 的 可 能 性 。 


居 移 动能 力 的 增长 缓慢 已 经 造成 了 非常 明显 
新 兴 硬 件 结构 之 间 发 展 的 不 平衡 和 不 
-个 “高 效 "的 数值 算法 应 该 使 单位 访 存 所 完成 的 浮 点 运算 次 数 
超 强 浮 点 运算 能 力 。 这 极 可 能 ; 
其 至 原则 上 的 根本 性 转变 。 本 文 尝试 面向 新 兴 计 算 机 体系 结构 ， 从 充分 释放 众 核 处 理 器 “ 匈 余 ” 计 
度 ) 的 广义 有 限 元 方法 ， 结 合 无 网 格 /粒子 /广义 有 限 元 
的 “计算 密集 度 可 调 * 的 共性 特征 ， 探 讨 硬件 效 能 和 数 


4b 


能 影 


ba 


响 


的 “ 浮 点 效率 鸿 
匹配 。 在 目前 的 新 兴 众 核 
尽量 加 大 ， 从 而 尽 可 能 
各 导致 数值 算法 设计 在 ， 


^I ^J 


也 享受 


s 
的 


E 


关键 词 : 百 亿 亿 级 计算 浮 点 效率 鸿沟 QA 协同 设计 广义 有 限 元 无 网 格 粒子 法 
1 百 亿 亿 次 计算 的 挑战 和 超级 计算 机 的 “ 浮 点 效率 鸿沟 ” 
表 1 目前 的 千 万 亿 次 系统 与 未 来 百 亿 亿 次 系统 〈 原 型 设计 ) 的 比较 . 
比较 内 容 泰坦 (Titan) 天 河 -II 百 亿 亿 次 系统 
(2011) (2013) (美国 DOE 方案 2020-2022 部 署 ) 
峰值 27Pflop/s 54.9Pflop/s 1.2Eflop/s 
ee 8.3MW (2Gflop/W) LT Sy ~20MW (50Gflops/W) 
i (1.935Gflop/W) 
710TB 1.4PB 
系统 内 存 ((82GB+6GB)x18688) ((64GB+3x8GB) 32-64PB 
x16000) 
3 " 1452Gflop/s 3431Gflop/s 
占星 
Hose (14141311) (42243009) Teron ona 
节点 内 存 232GB/s 304GB/s(?) 
带宽 (52+180) (64(?)+240) geet 
节点 16CPU+ 24CPU+ 
并 发 度 2688 Nvidia cores 3x57Xeon Phi cores O(10K) or O(1K) 
节点 总 网 6.35GB/s 
络 带宽 8GB/s (MPI broadcast) 200-400GB/s 
系统 节点 数 18688 16000 O(100K) or O(1M) 
系统 总 并 560640 3120000 
发 数 (299008 AMD cores+ (384000 Ivy Bridge cores O(1B) 
261632 Nvidia cores) + 2736000 Xeon Phi cores) 
平均 无 故障 
? ? 
时 间 ? ? O(«1day) 


民 据 2013 4F 11 


为 每 秒 107. (FIL) 次 浮 点 运算 )。 


A 一 


省 一 - 


月 的 超级 计算 机 500 3% (Top500) 排名 ， 中 国 的 天 河 开 号 (CPU-MIC 
异 构 架构 ) 成 为 日 前 世界 上 最 快 的 超级 计算 机 。Linpack 实测 
名 为 安装 在 美 


速度 为 54.9 petaflop/s(1 petaflops 
国 橡 树 岭 国家 实验 室 的 泰坦 (Titan， 


CPU-GPU 异 构架 构 )，Linpack 实测 速度 为 17.59 petaflop/s。 目 前 ,，“ 百 亿 亿 次 ”一 一 比 千 万 


亿 次 快 1000 倍 的 


超级 计算 技术 正在 成 为 世界 超 算 大 国 
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共同 挑战 的 目标 。 


超级 计算 面临 的 挑战 及 其 对 未 来 数值 算法 设计 的 可 能 影响 


美国 能 源 部 按照 电力 消耗 上 限 (如 不 应 超过 整个 胡 佛 水 坝 的 供电 能 力 ), 以 及 日 常 运 维 、 
果 有 费用 与 系统 造价 的 比例 关系 (如 保有 费用 不 应 超过 建造 成 本 )， 确 定 出 未 来 百 亿 亿 次 计 
算 机 设计 目标 , 其 中 总 体 功 耗 不 应 超过 20MWHl。 表 1 列 出 了 目前 世界 排名 第 一 和 第 二 的 天 
河 开 号 和 泰坦 (Titam) 以 及 美国 能 源 部 的 百 亿 亿 次 系统 原型 设计 的 主要 指标 。 通 过 表 1 中 的 比 
较 ， 我 们 可 以 发 现 两 个 事实 : (D) 天 河 开 号 的 峰值 为 泰坦 的 2 倍 ， 功 耗 也 正好 为 泰坦 的 2 倍 
多 一 点 ; (2) 两 台 机 器 的 功 耗 已 分 别 达到 8.3MW 和 17.8MW， 而 峰值 仅 分 别 为 百 亿 亿 次 的 
2.7% 和 5.49%。 事 实 (1) 说 明 : 以 天 河 II 号 和 泰坦 为 代表 的 目前 的 技术 进步 模式 是 线性 增 量 
式 的 ， 性 能 与 能 源 消耗 成 比例 提高 。 事 实 (2) 说 明 : 在 总 体 20MW 的 功 耗 约束 下 ， 目 前 的 线 
性 增 量 式 的 技术 进步 将 无 法 实现 百 亿 亿 次 计算 目标 ， 必 须 寄 希望 于 计算 技术 的 “颠覆 性 ” 进 
步 。 功 耗 效 率 成 为 了 实现 百 亿 亿 级 计算 所 面临 的 真正 挑战 1。 

另 一 方面 ， 在 大 规模 科学 与 工程 计算 
领域 ， 人 们 熟知 : 一 个 并 行 有 限 元 /有 限 差 ” 浮 点 效率 
分 程序 ， 即 使 可 获得 接近 线性 的 良好 并 行 Linpack 


加 速 ， 程 序 所 能 用 到 硬件 的 浮 点 性 能 也 只 - " 
3 x : 科学 及 工程 
占 到 峰值 性 能 很 少 的 一 部 分 。 无 论 是 以 单 应 用 
元 计算 为 主 的 显 格式 计算 〈 即 无 需 解 总 体 qm 
线性 方程 组 ， 如 差分 模板 〈stencil) 计算 )， N er 
还 是 以 线性 方程 组 求解 (核心 算法 多 为 稀 
AERE- t FE (Sparse Matrix-Vector multi- + 计算 访 存 比 
ply; spMV)) 为 主 的 隐 格 式 计算 ， 这 个 比 1 | 


例 一 般 很 少 超过 30%。 表 2 中 ， 我 们 针对 

差分 模板 计算 和 稀疏 算 阵 -向 量 乘 运算 ， 对 E. RT SEL SS SISTA” 
AMD Opteron 6274 fil Intel Ivy Bridge 两 球 最 新 多 核 CPU, 以 及 英 伟 达 k20x GPU 和 Intel MIC 
两 球 主流 众 核 处 理 器 的 “理论 浮 点 效率 上 限 ” 进 行 了 估算 。 不 考虑 时 间 局 部 性 ， 即 使 是 理论 上 
的 浮 点 效率 (假设 计算 数据 可 以 全 部 放 入 缓存 )， 结 果 仍 然 均 在 31% 以 下 。 另 外 ， 根 据 最 新 
Top500 排名 以 HPCG 《〈 和 采用 预 条件 共 斩 梯 度 法 求解 线性 方程 组 的 有 限 差分 计算 ) 对 现 有 


表 2 ”差分 模板 计算 和 迭代 法 求解 线性 方程 组 计算 的 理论 浮 点 效率 上 限 . 


AMD Intel nVidia Intel 
Opteron 6274 Ivy Bridge k20x GPU MIC 
峰值 (Gflop/s) 141 422 1311 3009 
峰值 带宽 比 (m) 21.69 52.75 58.26 100.3 
3 点 差分 (3D) 6.5% 2.7% 2.4% 1.4% 
5 点 差分 (3D) 14.4% 5.9% 5.3% 3.1% 
7 点 差分 (3D) 22.3% 9.2% 8.3% 4.8% 
9 点 差分 (3D) 30.2% 12.4% 11.2% 6.5% 
spMV 9.22% 3.79% 3.43% 1.99% 
超级 计算 机 的 测试 结果 表明 ，( 即 使 相对 LINPACK 实测 值 ) 最 好 的 浮 点 效率 仅 为 4.1% (日 本 
的 “ 京 ” 超 级 计算 机 ) 〈 表 37。 然而 ， 用 测试 标准 程序 Linpack 测 得 的 超级 计算 机 的 浮 点 效率 


一 般 可 以 达到 80% 以 上 。 日 本 的 “ 京 "超级 计算 机 的 Linpack 测试 结果 甚至 高 达 93%。 实 际 科 
学 与 工程 应 用 的 浮 点 效率 与 超级 计算 机 的 Linpack“ 名 义 ” 浮 点 效率 之 间 形 成 了 一 个 明显 的 
“ 浮 点 效率 鸿沟 ”( 表 3)。 


1 从 长 远 看 ， 成 功 迈 向 百 亿 亿 级 计算 的 真正 意义 ， 绝 不 仅仅 是 完成 一 台 百 亿 亿 次 计算 机 的 建造 ， 它 将 触发 整个 计算 技术 产业 ， 包 括 从 超级 计算 机 到 
手持 设备 的 颠覆 性 技术 进步 。 比 如 ， 百 亿 亿 级 时 代 ， 你 的 手机 处 理 能 力 将 更 强 、 待 机 时 间 却 更 长 。 

? 加速 比 为 一 个 并 行程 序 /并 行 算法 的 并 行 性 能 的 主要 评价 指标 。 在 问题 规模 一 定 的 情况 下 ， 使 用 两 个 处 理 器 核 应 该 获得 比 一 个 处 理 器 核 快 一 倍 的 计 
算 速 度 ， 即 理论 加 速 比 应 为 2， 以 此 类 推 ,， 对 n 个 处 理 核 ， 理 论 加 速 比 应 线性 增加 到 n。 一 个 可 获得 线性 加 速 的 并 行程 序 / 算 法 即 被 认为 是 很 好 的 并 
行程 序 /并 行 算 法 。 
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因此 我 们 目前 面 对 的 现状 是 : 一 方面 制造 百 亿 亿 次 超级 计算 机 面临 越 来 越 大 的 能 耗 挑 
馈 级 计算 机 的 实际 浮 点 效率 却 


表 3 现 有 超级 计算 机 在 真实 应 用 下 的 浮 点 计算 效率 (数据 来 源 : Jack Dongarra) 
x Linpack "m 
ids SERO 计算 核 su Top500 真实 应 用 峰值 真实 应 用 峰值 
数目 排名 (pflops) /Linpack 峰值 
(pflops) 
中 国 广 州 天 河 -2 
Xeonl2c22GHz 3,120,000 33.9 1 0.580 1.7% 
超 算 中 心 Intel Xron Phi 57c+ 定 制 核 
re AC 超级 计算 机 
Bt 超级 计 
理化 高 级 计 “ 窗 二 通 SPARCG4VIlltxsc 705,024 — 10.5 4 0.427 4.1% 
算 研究 院 ee 
美国 能 源 部 泰坦 
橡树 岭 国家 Nw Kepiercpu tac 540,640 — 17.6 2 0.322 1.8% 
实验 室 + 定制 核 
美国 能 源 部 
Mira BlueGene/Q 
阿 贡 国家 实 。 Power pcc 786,432 8.59 5 0.101 1.2% 
FE Hil AZ 
瑞士 国家 超 Piz Daint 
Cray XC30,Xeon 8C 115,984 6.27 6 0.099 1.6% 


级 计算 中 心 


产生 “ 浮 点 效率 鸿沟 ”的 原因 
目前 大 型 线性 方程 组 求解 主要 采 月 
心 多 归结 为 稀疏 外 


flop/word, 


= 


标 


+Nvida Kepler 14C+ 定 制 核 


E 阵 -向 量 乘 ，T 
word 为 浮 点 字 长 ， 为 4 或 8 字 节 )， 是 典型 的 访 存 密集 型 运算 (通常 人 们 笼统 地 
称 科学 与 工程 计算 为 计算 密集 型 ， 是 与 其 它 的 计算 机 应 月 
八 测 试 所 对 应 的 算法 是 稠密 矩阵 的 直接 法 求解 ， 算 


并 不 复杂 : 科学 与 工程 计算 相当 一 
Hoe (Krylov) 子 空间 进 代 法 "， 这 类 算法 最 终 的 核 


部 分 为 求解 线性 方程 组 。 


条 稀 芷 矩阵 -向 量 乘 算法 的 “计算 访 存 比 "不 超过 2: 1 CHUTE: 


precision GEneral Matrix Multiply (DGEMM)), “H 
届 于 计算 密集 型 运算 。 WA, PX 


EE 


而 且 坏 消息 是 ， 这 种 差距 目 


种 本 质 性 的 差异 直接 决定 了 Linpack 测 试 和 实际 科学 与 了 


类 型 相 比 较 而 言 的 )。 而 Linpack 
法 核心 为 双 精 度 和 矩阵 相 乘 (Double 
EWT” H n: 1 (Cn ARBRE), 


前 还 在 进一步 扩大 。 


2 ” 众 核 计算 一 一 进一步 扩大 的 “ 浮 点 效率 鸿沟 ” 


2 数据 源 
样 的 数据 , 我 们 对 
的 双 精 度 浮 点 峰 


C2050 哪 


图 2 中 左 图 为 7 


款 产品 更 好 ? ”PC 时 


RIAA A, j 


Tesla C1060 和 Tesla C2050 y 


EK GPU 在 进 


法 的 “计算 访 存 比 " 具 有 本 质 性 差别 。 而 正 是 算法 上 的 这 
[ 程 应 用 之 间 在 浮 点 效率 上 的 明显 差 


HT on AXE GPU 相对 于 多 核 CPU 的 性 能 优势 。 基 于 同 


行 “ 双 精度 稠密 矩阵 -矩阵 乘 "运算 时 的 性 和 


现款 众 核 GPU 处 理 器 做 一 个 横向 比较 。C2050 
让 速度 是 C1060 的 6.6 倍 。 如 果 我 们 问 :“ 对 于 科学 与 工程 计算 ，C1060 和 
尺 的 经 验 是 峰值 速度 越 高 越 好 。 


EE 差异 。 计 算 通 过 


调用 cuBLAS3.1 库 中 的 DGEMM 子 函 数 实现 。DGEMM 正 是 超级 计算 机 性 能 评测 程序 
应 用 程序 和 一 个 稠密 矩阵 计 


Linpack 的 核心 算法 。 图 2 中 右 图 则 是 针对 13 个 真实 科学 与 工科 
其 中 大 多 数 为 有 限 元 /有 限 差分 计算 一 一 的 测试 结果 。 这 些 真 实 应 用 的 核心 算法 


算 程序 


DI fii HE MAE [6] HEE RAT] 


WERKA KIRE EIA CG， 广 义 最 小 残 量 法 GMRES, PAHAT REV BiCG 等 。 
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r1 


看 到 ， 对 于 Linpack 测试 而 言 ，C2050 速度 是 C1060 的 两 倍 ， 


而 对 于 真实 应 用 ，C2050 仅仅 比 C1060 快 了 约 4096 


超级 计算 面临 的 挑战 及 其 对 未 来 数值 算法 设计 的 可 能 影响 


对 于 Linpack 测试 ，C2050 和 C1060 


Gflops MISERERE Gflops 
200 2d 双 精 度 计 算 
—H- Tesla C2050(ECC 关 ) 
160 16 -€-- Tesla C2050(ECC 开 
-m- Tesla C2050(ECC 关 ) -0- Tesla C1060 
-@- Tesla C2050(ECC3F) —9- Intel Xeon 5550 
120 -O- Tesla C1060 12 
—— Intel MKL 4 线程 
80 8 
40 4|8 
Oe xn x 9? 9 gm 0 aana AEREA ENE 
Q 5 a Oo 3 2 iom m8 LOS oL oL oL on o4 EE ROS 
x x — CN st co St a4 路 o BK iW mE tb ip R & 
Oo N x x x x R ® z m z Eg oa R 
LO 一 十 co CO N ES ig K HB qo Eg 中 
Nn ow A + Oo OD Ip FP BR ti 
Oo ò à = ES ic X g 
— N sF co 十 ur 
矩阵 大 小 
(a) 双 精 度 稠密 和 矩阵 乘法 DGEMM (b) 双 精 度 稀 疏 矩 阵 -向 量 乘 
图 2，Tesla C2050 和 Tesla C1060 两 款 众 核 GPU 处 理 器 的 性 能 测试 


原始 数据 来 源 : 英 伟 达 公司 网 站 外。 对 于 Linpack 测试 : C2050 re C1060 性 能 提升 显著 ; 
对 于 实际 的 科学 与 工程 计算 : 前 者 无 特别 明显 的 “能 效 ” 优 势 。 


性 能 差异 明显 


， 而 对 于 真实 应 用 ， 二 者 性 能 则 


X4 C2050 和 C1060 GPU 的 主要 性 能 指标 


ab 
非常 接近 。 GPU 功 耗 峰值 内 存 带 宽 
造成 这 种 性 能 表现 上 的 差异 的 厦 因 如 下 ; C1060  188W  78Gflops 102GB/s 

成 这 种 性 能 表现 上 的 差异 的 原因 如 下 C2050 247W 515Gflop/s 144GB/s 


DGEMM 是 计算 密集 型 操作 , 其 性 能 提升 理论 


上 由 人 硬件 的 浮 


也 受到 GPU 寄存 器 大 小 的 限制 ， 使 得 分 块 矩 阵 n 120 
不 能 很 大 ， 大 大 降低 了 实际 计算 的 计算 访 存 比 )。 。 100 


点 计算 能 力 的 提升 决定 (但 实际 上 ， 即 使 这 一 计算 密集 型 算法 ， 在 C2050 上 


而 稀疏 矩阵 -向 量 乘 是 典型 的 访 存 密集 型 操作 ， 划 
性 能 提升 主要 由 硬件 的 访 存 带宽 的 提升 幅度 决 
定 。C2050 的 访 存 带宽 (144GB/s) Lk C1060 60 
(102GB/s) 高 出 46%。 这 正好 解释 了 为 什么 对 于 40 
真实 应 用 ，C2050 仅仅 比 C1060 HET 40%。 实 际 


80 


上 ， 如 果 考 虑 功 耗 、 性 价 比 ， 特 别 是 功 耗 因素 ， " 

C1060 5j C2050 相 比 ， 对 于 科学 与 工程 计算 而 言 ET 

可 能 并 不 算 差 。 如 果 从 浮 点 效率 〈 能 耗 效率 ) 来 9 3 9 = 

fi, C1060 (14%) 反而 比 C2050 (3%) 更 高 5 2 $ E 

虽然 C2050 比 C1060 浮 点 性 能 快 6.6 倍 ， 但 效率 ô £ B 

相 比 C1060 却 大 大 降低 。 3 ^ 

图 3 进一步 综合 性 地 给 出 了 典型 应 用 算法 在 目前 图 3， 处 理 器 性 能 与 典型 应 用 

四 款 主 流 处 理 器 上 的 浮 点 效率 的 变化 趋势 。 我 们 的 浮 点 效率 变化 曲线 
发 现 ， 由 于 真实 科学 与 工程 计算 的 访 存 密集 的 特 由 于 真实 科学 与 工程 计算 的 访 存 密 
征 所 限 ， 它 们 并 没有 简单 地 从 浮 点 性 能 的 单一 增 、 集 的 特征 所 限 ， 它 们 并 没有 简单 地 从 浮 


长 中 获 益 ， 反 而 出 现 浮 点 速度 越 快 ， 浮 点 效率 越 


低 的 怪 现 象 。 


浮 点 效率 低下 即 计 算 能 力 的 浪费 。 在 “ 功 耗 "成 为 百 亿 亿 次 


点 性 能 的 单一 增长 中 获 益 ， 反 而 出 现 浮 


点 速度 越 快 ， 浮 点 效率 越 低 的 怪 现 象 。 
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HER, AINE AR BTA)” RAT BUR R 


RIK 


E 要 的 现实 意义 。 


3 ”难点 所 在 ;传统 数值 算法 与 新 兴 众 核 硬 件 结构 之 间 的 不 匹配 


“ 浮 点 效率 鸿沟 ”的 存在 , 不 是 计算 机 硬件 设计 不 
好 ,在 过 去 40 年 中 , CPU 核心 的 速度 提高 了 近 千 倍 ， 
但 是 内 存 的 读 写 速度 提高 却 不 大 ; 更 为 关键 的 是 ， 
内 存 的 延迟 没有 降低 。 如 从 DDR 到 DDR3， 内 存 的 
传输 速度 成 倍 地 提高 ， 但 是 内 存 延迟 没有 减少 ， 反 


而 还 有 所 增加 。 这 使 得 内 存 延迟 已 经 成 为 了 制约 


CPU 实际 处 理 能 力 提高 的 瓶颈 。 在 实际 应 用 


H, CPU 
始终 都 在 等 待 内 存 访 问 。 当 内 存 端口 满 负 载 时 ，CPU 


核心 有 超过 $0% 时 间 是 空闲 的 。 这 就 是 早 在 1994 年 
便 被 提出 的 “内 存 墙 (Memory WalD)”DI 问 题 。 目 前 ， 非 


易 失 性 存储 器 CNVM) 以 及 三 维 封 装 技术 〈3D 
Integration) 可 互 成 为 具有 “颠覆 性 ”的 新 一 代 存 储 技 
术 。 但 是 无 论 采 取 何 种 快速 存储 技术 , 计算 机 的 冯 : 诺 


根本 的 原因 。 克 服 内 存 墙 将 是 一 个 长 期 、 


计算 密集 度 O(1) 


伊 曼 “存储 程序 原理 ”架构 可 能 是 造成 内 存 墙 问题 最 
巨大 的 技术 挑战 由。 


稠密 


pteron 6274 FR 


O(log(N)) 
—_ L5 


Vol.12 No.4 
Oct. 2014 


Intel Ivy Bridge 
nVidia k20x GPU 
Intel MIC 


| 


=> AMDO 


CC 


图 4， 性 


与 访 存 带 宽 的 发 展 趋势 


线性 代数 /LINPACK(BLAS3) 


计算 访 存 比 TOC ale 


图 5， 处 理 器 结构 与 数值 算法 在 基本 计算 特征 上 
“alg” 为 浮 点 格式 字 长 ， 单 精度 为 4， 双 精度 为 8 


“内 存 墙 ?问题 和 新 兴 多 核 / 众 核 技术 的 发 展 正在 使 得 处 至 


) 


硬件 的 计算 性 能 、 数 据 移 动 性 能 之 比值 (计算 机 学 科 ) 
CE 
众 核 /0((5-10)xalg) 


v= 
BIOS) — Er 


的 不 匹配 和 不 协调 


器 计算 能 力 和 访 存 能 力 之 间 的 


差距 愈加 扩大 喇 。 如 图 4 所 示 , 我 们 同样 针对 AMD Opteron6274 和 Intel Ivy Bridge 两 款 最 新 


多 核 CPU， 以 及 英 伟 达 k20x GPU 和 Intel MIC 两 


增长 趋势 。 从 中 可 以 看 出 ， 处 理 器 性 能 呈现 3 


^ 内 存 延迟 是 指 内 存在 接 到 访问 命令 后 ， 要 等 一 段 时 间 ， 才 能 伟 


ay, 
IA 


输 数 据 的 速度 ) 有 所 提高 ， 但 是 DRAM 的 延迟 不 但 没有 什么 改善 ， 


反 
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最 新 众 核 处 到 


到 并 


可 数据 。 这 个 延迟 超过 100 纳 秒 。 过 去 30 年 ， 
而 有 继续 扩大 的 趋势 ”。 


器 ， 画 出 峰值 与 访 存 带宽 
指数 级 增长 ， 而 访 存 能 力 几乎 呈 线 性 增长 或 微弱 


DRAM 内 存 的 传输 速度 〈 仅 仅 是 传 


超级 计算 面临 的 挑战 及 其 对 未 来 数值 算法 设计 的 可 能 影响 
的 指数 级 增长 ， 二 者 间 的 差距 也 呈 指 数 级 增长 。 实 际 上 ， 从 单 核 、 多 核 到 众 核 , 处 理 器 的 “ 计 
算 访 存 能 力 比 值 " 越 来 越 大 ,已 经 发 生 了 约 一 个 数量 级 的 改变 ( 见 图 5 中 计算 访 存 比 的 标注 )。 


“计算 访 存 能 


力 比 值 ”高 意味 着 处 理 器 更 加 喜欢 高 “计算 密集 度 ” 的 数值 算法 ， 即 和 


位 访 


存 所 进行 的 浮 点 运算 次 数 很 多 的 数值 算法 。 然而 , 数值 算法 或 应 用 数学 学 科 的 发 展 却 是 非常 


缓慢 的 一 一 与 计算 机 处 型 


数值 方法 的 计算 密集 度 进行 归 类 ， 我 们 会 发 现 : 
仍然 以 “线性 计算 密集 度 ” 算 法 为 多 ， 密集 地 集 ! 


器 技术 激进 演化 的 发 展 节奏 完全 不 同 。 从 


图 5 可 以 看 出 ,如 果 按 照 
目前 80-90% 的 科学 与 工程 计算 的 核心 算法 


却 远 远 位 于 第 头 的 右 端 ,处 于 


现 明显 的 不 适应 和 不 协调 。 


4 ”研究 现状 及 发 展 趋势 一 


综合 上 述 表 1、 


图 3 和 图 


4， 我 


们 可 以 归纳 出 以 下 基本 事实 : 


(1) 


= 


型 真实 应 用 的 浮 点 效率 远 
远 小 于 Linpack XK OFF 


DA f iit kB. BE- p] E 36 AJ FA À 


算法 的 线性 
率 均 低 于 10%) 


» M 


方程 组 求解 ， 效 
TEL, 峰 


值 


带宽 比 越 高 ， 实 际 
点 效率 反而 越 低 。 


(2) 


组 克 雷 洛 夫 子 空 
解 的 核心 计算 ) 
EG? il a, 为 最 典 


(3) 理论 


到 改善 。 


REM- ER CRE 
HARIR 
“计算 访 存 
型 的 访 存 密 
集 型 计算 ， 浮 点 效率 最 低 。 


上 ， 差 分 模板 计算 的 效 
率 可 以 通过 改变 差分 格式 得 


事实 (1) 说 明 ， 处 理 器 的 峰值 


应 用 的 浮 


方程 


或 峰 


于 箭头 的 左 端 。 而 便 件 的 计算 访 存 能 力 比 值 


里 器 技术 的 激进 快速 发 展 与 相对 发 展 缓慢 的 数值 算法 之 间 已 经 出 


“好 ”的 数值 算法 的 定义 在 发 生 改变 


通过 提高 数值 方法 的 计算 密集 度 ， 在 充分 发 挥 众 
应 用 问题 数值 求解 精度 


协同 设计 区 间 : 
核 “ 元 余 ” 的 浮 点 运算 能 力 的 同时 ， 提 高 


iF 7 REL 

数 信 粮 启 | 学 与 工程 计算 i 
Linpack 导 向 下 

90% 的 硬件 设计 


80% 


a" | 缩短 效率 鸿 党 
20% 

0 集 度 

PC/ 单 核 众 核 

科学 与 工程 计算 ”科学 与 工程 计算 LINPACK 

EIS ae areata 稠密 矩阵 

低 阶 方法 高 阶 方法 全 局 支持 域 方法 
(RBF， 谱 方法 等 ) 


Á 广义 有 限 元 、 粒 子 方法 等 
(PCB BRTRIEE? 
\ 设 计 惯性 思维 方向 \ 


构 的 协同 设计 方向 


值 带宽 比 越 高 ， 实 际 应 用 的 浮 点 效率 N 

越 低 ， 反 映 了 正在 加 剧 的 “ 浮 点 效率 图 6， 面 向 新 兴 处 理 器 体系 结构 

鸿沟 ”问题 ， 反 映 了 计算 机 硬件 与 数 的 数值 算法 的 协同 设计 

值 算法 之 间 的 不 适应 和 不 协调 。 事 实 (2) 表 明 ， 需 要 针对 新 兴 计算 机 体系 结构 ， 在 算法 层面 
对 线性 方程 组 的 求解 方法 进行 改变 或 重新 选 优 ， 以 避免 过 低 的 浮 点 效率 一 在 目前 的 多 核 / 
众 核 处 理 器 上 ， 如 果 关 键 数 值 算法 不 作出 改变 , 浮 点 效率 将 无 法 得 到 真正 有 效 的 改善 。 事实 


(3) 揭 示 了 克服 面向 浮 点 效率 问题 


时 ， 我 们 可 以 采取 的 


\ 体 措施 一 一 改变 数值 算法 。 通 过 改 


变 差分 模板 的 阶 次 ( 低 阶 党 高 阶 


)， 在 差分 阶 次 或 差分 精度 提高 的 同时 ， 浮 点 效率 也 得 到 了 改 


善 ( 计 算 效 率 和 算法 精度 的 双赢 )。 这 说 明 “ 好 ”算法 的 定义 正在 发 生 改 变 。 


过 去 被 认为 是 计算 “费力 ”的 高 精度 算法 可 
百 亿 亿 次 系统 上 , 科学 与 工程 应 用 算法 应 该 使 单位 访 存 所 完成 的 浮 点 运算 次 数 最 大 化 ， 从 而 


地 说 ,在 未 来 


能 是 未 来 的 “好 ”算法 。 或 者 简单 
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尽 可 能 地 享受 新 的 “免费 午餐 ”一 一 浮 点 计算 能 力 。 在 数值 方法 设计 方面 ， 在 PC 时 代 ， 


能 
浮 点 运算 能 力 受 限 ， 算 法 设计 总 是 喜欢 “线性 计算 密集 度 ” 的 算法 ， 喜 欢 稀 玻 性 。 如 图 6 


所 示 ， 过 去 我 们 沿 “ 向 左 ” 箭 头 方向 思考 。 而 面 对 众 核 的 高 浮 点 运算 能 力 和 数据 移动 能 力 的 相 


对 不 足 ， 未 来 应 该 强调 矩阵 “适当 ”稠密 ， 增 加 计算 密集 度 ， 鼓 励 沿 “ 向 右 ” 箭 头 方向 思考 。 
将 带 来 科学 与 工程 应 用 算法 设计 在 思路 上 甚至 原则 上 的 根本 性 转变 。 


5 ”面向 众 核 计算 的 数值 方法 协同 设计 


这 


按照 图 6 的 预测 : 无 网 格 、 粒 子 类 方法 以 及 (单位 分 解 ) 广义 有 限 元 等 一 类 过 去 被 认为 
计算 “费力 ”的 但 具有 更 好 的 数据 局 部 性 的 新 兴 数 值 方法 值得 我 们 从 发 挥 众 核 计算 能 力 的 角 


度 重 新 给 予 关注 。 


传统 上 ， 数 值 方法 研究 关注 的 只 是 精度 、 收 敛 性 等 数学 
指标 ， 从 来 都 不 会 去 关注 “计算 密集 度 ” 这 一 来 自 于 计算 机 新 
兴 体 系 结构 对 数值 算法 的 新 要 求 。 受到 表 2 中 有 限 差 分 浮 点 
效率 变化 的 启发 ,图 7 对 无 网 格 /粒子 类 方法 进行 了 类 似 的 考 
查 。 假 设 胞 元 cell) 中 的 平均 粒子 数 为 n。 中 心 胞 元 内 粒子 
的 重用 近似 等 于 (或 小 于 ) 周围 胞 元 内 粒子 总 和 ,为 GB, 
d 为 问题 维 数 。 周 围 胞 元 内 粒子 的 重用 大 臻 等 于 (或 小 于 ) 


o Ogle C 0O 9o 
9 © 090909Q 9 © 
9o 9ojlo © 9o Og 
IC o 0j9o OL o0 
o 9 gjo oO 9o 9o 
IQ o W019 o Of Qo O0 


EK nB) /34。 


XS 无 网 格 /粒子 类 方法 “单元 ”计算 的 理论 浮 点 效率 
Cell 中 粒子 数 AMD Opteron 6274 Intel Ivy Bridge nVidia k20x GPU Intel MIC 
2x2x2(n=8) 71.0% 29.2% 26.4% 15.3% 
3x3x3(n=27) 100.0% 98.5% 89.2% 51.8% 


心 胞 元 内 粒子 个 数 ， 为 n。 这 一 类 方法 的 “平均 ”计算 密集 ”图 7， 无 网 格 /粒子 类 方法 


对 这 一 类 数值 算法 ， 我 们 可 以 针对 前 面 列 出 的 四 款 处 理 器 估算 出 理论 浮 点 效率 。 如 表 5 
所 示 ， 我 们 发 现 ,在 同样 处 理 器 上 ， 通 过 提高 算法 精度 ， 浮 点 效率 也 大 幅 提 升 。 理 论 分 析 表 
明 ， 这 类 无 网 格 /粒子 类 方法 以 及 前 面 提 及 的 有 限 差 分 方法 ， 一 个 最 大 的 共性 特点 是 ， 通 过 


改变 算法 格式 (差分 阶 或 胞 元 /影响 半径 的 大 小 〉 可 以 主动 控制 处 理 器 的 浮 点 效率 。 


受 此 启发 ， 我 们 基于 前 期 发 展 的 无 额外 自由 度 的 广义 有 限 元 方法 5 下 工 9， 探 讨 在 众 核 


计算 平台 上 计算 密集 度 调 节 的 重要 性 和 意义 。 下 面 先 简单 介绍 这 一 新 的 广义 有 限 元 方法 。 


无 额外 自由 度 的 广义 有 限 元 方法 
单位 分 解 〈partition of unity) 广义 有 限 元 方法 的 研究 一 般 认为 源 于 梅 伦 克 (Melenk) 


和 


巴 布 什 卡 〈Babugka) 的 单位 分 解 方法 0 和 单位 分 解 有 限 元 方法 C〈 类 似 的 思想 又 可 追溯 到 
该 作者 早期 的 工作 P2)。 几 乎 同时 ， 杜 阿尔 特 (Duarte) AISLE (Oden) PPAP aR LAG DT 


者 石 根 华 提出 的 数值 流 形 方法 JJ。 单位 分 解 函 数 概 念 本 身 非常 简单 ， 既 可 使 用 无 网 格 扣 


法 ， 称 为 hp-cloud 方法 。 在 思想 上 非常 相近 但 发 表 时 间 更 早 的 工作 ， 可 见于 由 我 国旅 美学 


ER 


方法 构造 ， 也 可 直接 使 用 有 限 元 的 形 函 数 。 作 为 有 限 元 方法 的 自然 延伸 ， 基 于 有 限 元 形 函 数 
的 单位 分 解 方法 获得 了 极 大 的 关注 和 成 功 的 应 用 ， 如 斯 特 劳 鲍 里 斯 (Strouboulis) 和 巴 布 什 
卡 89B7 和 杜 阿尔 特等 深入 发 展 的 广义 有 限 元 "3%。 为 了 清晰 指 代 ， 下 面 的 叙述 中 “广义 有 


限 元 法 (GFEM)” 主 要 指 代 基于 有 限 元 形 函 数 的 单位 分 解 方法 。 


> Numerical Manifold Method, NMM 


超级 计算 面临 的 挑战 及 其 对 未 来 数值 算法 设计 的 可 能 影响 


广义 有 限 元 的 核心 是 单位 分 解 逼 近 : 
u" (x) >) Nu + i Ni Aliw (1) 


tB Y N(x) = 1 FRET 4k Ef EAM gi 为 节点 i 支撑 域 上 的 用 户 自 定义 局 部 近 
似 函 数 ，ai)，aio)，.… 为 局 部 近似 函数 引入 的 节点 i 的 广义 自由 度 或 额外 自由 度 。 


在 现 有 的 广义 有 限 元 方法 中 ， 节 点 上 自由 度 个 数 会 随 独 局 部 近似 函数 的 阶 次 变化 而 变 
化 ， 因 此 提高 局 部 函数 的 阶 次 不 会 改变 该 方法 的 计算 密集 度 〈 见 稍 后 图 12 的 说 明 )。 下 面 ， 
我 们 提出 一 种 “无 广义 自由 度 的 广义 有 限 元 方法 ”( 下 称 “ 新 方法 ”)。 


假设 PP 为 环绕 节点 i 的 所 有 单 
元 组 成 的 “单元 片 (nodal patch) ”, 
r 代表 单元 片 尺 寸 。 根据 网 格 特征 ， 
单元 片 尺寸 可 通过 两 种 方式 定义 : 
对 于 规则 结构 化 网 格 或 非 均 匀 网 
格 ， 为 环绕 节点 i 的 mel 层 单元 的 


文 撑 域 的 大 小 ( 见 图 8(a)); 对 于 任 (a) 结构 化 网 格 (b) 任意 网 格 
意 均匀 网 格 ， 则 可 简单 取 为 节点 i 图 8， 单 元 片 定义 


的 影响 圆 (或 球 ) 半径 ( 见 图 8(b))。 
节点 宇 定 义 为 “ 片 星 Cpatch star) ”。 若 无 特别 说 明 ， 指标 i 做 特意 保留 ， 特 指 “ 片 星 节 点 (patch 
star)”， 指 标 了 Gai) 表示 “ 非 片 星 节 点 (non-patch star) ”。 {x Ix, € Pi | 表示 单元 片 P 上 
所 有 节点 的 集合 。 


首先 ， 使 用 节点 Lx, [xy © pr] ORE un), x, e pr 


元 片上 的 一 个 局 部 近似 : 


Fr 


u (x) - et (x), Q) 
gi (x) SEC, mi NEM Er AC, an 为 节点 未 知 数 。 
将 局 部 近似 wx 直接 取代 标准 有 限 元 的 节点 位 移 us 


N 
4| 


w^ (x) Ni (x) (3) 
则 可 得 出 一 种 新 的 逼近 形式 
"Bw (Seem (4) 


插值 (4) 与 插值 (3) 均 基于 同一 有 限 元 网 格 的 节点 构造 ， 并 没有 引入 新 的 节点 或 自由 度 。 
将 式 (4) 中 的 项 展开 、 然 后 按 相 同 节 点 归 类 、 重 新 组 合 后 可 得 到 
u" (x) 223 (Ng 1» sd Ng In (5) 
其 中 gr 为 在 单元 片 k 上 节点 i 的 局 部 函数 (注意 gr 和 gi EREM ni foo Pr EAA 
数 ， 或 包含 节点 i 的 所 有 单元 片 的 数量 ，N, 表示 新 的 形 函 数 。 


比较 (4) 和 (1)， 可 以 清楚 地 看 出 ，(4) 实 质 上 也 是 一 种 单位 分 解囊 近 。 然 而 ， 新 单位 分 解 
逼近 具有 两 个 独一无二 的 特点 : (a) 不 包含 任何 广义 自由 度 ; (b) 只 要 局 部 逼近 函数 在 各 自 的 
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“ 片 星 ” 处 插值 , 则 无 论 其 在 “ 非 片 星 ?处 播 值 与 否 , 最 终 构造 出 的 逼近 函数 一 定 是 全 局 插值 的 。 
新 方法 的 特性 逐条 详细 讨论 如 下 : 


(1) L 有 现 有 广义 有 限 元 方法 的 高 阶 性 质 。 通 过 增 大 单元 片 尺寸 ， 可 以 得 到 “高 阶 ” 的 局 
部 函数 ， 从 而 实现 高 阶 插值 。 


《2) 不 包含 广义 自由 度 ， 并 且 待 求 的 总 自由 度 不 会 随 着 局 部 函数 的 阶 次 改变 而 改变 。 这 
是 现 有 广义 有 限 元 方法 不 具有 的 特点 , 这 直接 导致 线性 无 关 性 Clinear independence? 
和 与 标准 有 限 元 一 样 的 稳定 性 5 。 这 两 点 特性 是 与 现 有 广义 有 限 元 方法 最 大 、 最 本 
质 的 区 别 。 


G) RAEE. 局 部 函数 可 以 采用 某 种 插值 方法 , 如 拉 格 朗 日 插值 或 径 向 基 函 数 等 ， 
或 者 采用 具有 ” 单 点 插值 ?特性 的 某 种 逼近 方法 来 构造 。 对 于 后 者 ， 可 以 使 用 移动 最 
小 二 乘法 来 构造 局 部 函数 ， 只 需 强 制 该 最 小 二 乘 拟 合 函数 在 “ 
该 “ 单 点 插值 "约束 很 容易 满足 。 在 这 些 情况 下 ， 所 构造 出 的 全 局 双 近 函数 自然 共有 
全 局 插值 性 质 。 


5.1 局 部 函数 的 构造 一 一 “ 单 点 插值 ?最 小 二 乘 逼近 


局 部 函数 的 构造 是 基于 单元 片上 的 节点 集 来 实现 的 。 现 有 的 标准 化 技术 ， 如 无 网 格 方 
法 5 并 3 或 者 网 格 类 方法 ， 都 可 借用 过 来 。 唯 一 的 要 求 是 : 如 果 使 用 逼近 格式 ， 比 如 移动 
最 小 二 乘 后， 它 应当 而 且 只 需 在 单元 片 的 “ 片 星 ”处 具有 插值 性 质 。 对 于 结构 化 网 格 ， 可 以 使 
] 拉 格 朗 日 多 项 式 插值 ， 对 于 任意 网 格 ， 可 以 使 用 径 向 要 函数 [的] 构造 局 部 函数 。 这 些 构 
造 过 程 直 观 人 简单 ， 这 里 不 再 次 述 。 


下 面 重 点 介绍 我 们 提出 的 适用 于 任意 网 格 的 “ 单 点 插值 ”移动 最 小 二 乘 逼 近 ! 9 。 这 一 
方法 使 我 们 能 够 使 用 局 部 逼近 来 构造 出 全 局 插值 的 广义 有 限 元 方法 。 我 们 仍然 用 指标 宇 特 指 
“HATE”, j 表 示 “ 非 片 星 节 点 *。 在 以 i 为 中 心 的 单元 片上 ， 局 部 双 近 定义 为 

u, (x)=p"(x)a(x),x € P; (6) 
FERAL p 可 以 包含 解 的 先 验 知识 ， 对 于 光滑 问题 ， 可 以 直接 取 为 任意 高 阶 多 项 式 ， 对 
具有 局 部 特征 的 问题 ， 如 裂纹 尖端 ， 可 以 是 非 多 项 式 特殊 函数 一 简 言 之 ， 基 函数 是 可 以 看 
据 竺 求 问题 解 的 先 验 知识 进行 自 定 义 的 , p 的 长 度 与 单元 片 尺寸 ch 相 适 应 (h 为 网 格 尺寸 )， 
p 也 可 以 包含 解 的 先 验 知识 ， 如 扩展 有 限 元 法 “中 使 用 的 裂纹 尖端 基本 解 等 ，a 为 待定 系数 。 

为 了 构造 一 个 在 “ 片 星 ”处 插值 的 移动 最 小 二 乘 逼 近 , 我 们 使 用 下 面 的 带 有 约束 条 件 的 最 
小 二 乘 逼近 


= 并 


其 


本 = 二 ALTE (x)(p'( x,)a - u,) +A(p"(x;)a—u; ) (7) 


其 中 严 是 单元 片 主 上 的 节点 总 数 。 为 了 便于 程序 实现 ,“ 片 星 节 点 完 也 包含 在 求 和 项 中 ， 这 
样式 (7) 的 第 一 部 分 与 标准 移动 最 小 二 乘 器 完全 一 样 ,， ?为 拉 格 朗 日 乘 子 , 用 于 强制 u(x) 精 确 
满足 “ 单 点 插值 "条件 


p'(x;)a -u, (8) 


n di E Ee e — FR LAL 


对 最 小 二 乘 范 数 取 极 值 可 得 到 “8 


c 


® Extended Finite Element Method , XFEM 


超级 计算 面临 的 挑战 及 其 对 未 来 数值 算法 设计 的 可 能 影响 


u (x)= > ø (x) uy (9) 

k=l 

I 
r E 1 dI 1 E 
h: (x) 三 p. (sa Pr -AT AoA Pe alba) (10) 
m 1 

A=) WPP: (11) 

要 
T p y - Ji : 12 
p E PME (12) 


Aq JATKE, ATA Ai 的 第 一 个 元 素 ，6 AGES ARAL (Kronecker delta). uf 
以 证 明 ， 当 x=x; 时 ， 将 p7 =[1,0,0,...] ATI g (x;)=1， 从 而 说 明 局 部 函数 在 “ 片 星 节 
点 ”处 满足 插值 性 质 。 

在 公式 (7) 中 , A we1,“ 单 点 插值 ”移动 最 小 二 乘 便 退 化 为 “ 单 点 插值 ”最 小 二 乘 (SILS)。 
相对 “ 单 点 播 值 移 动 最 小 二 乘 ,“ 单 点 插值 最 小 二 乘 (特别 是 其 导数 ) 具有 计算 速度 快 的 优 


yy o 


使 用 SIMLS / SILS fE29J5 SEE, 30) SCA Pe oci n] WA eae a TF 


ua) Som Sa | (13) 


IEP gg 为 单元 片 Pr 上 的 SIMLS 或 SILS 局 部 函数 。 


SIMLS/SILS 适用 于 任何 维 数 的 规则 或 不 规则 网 格 。 基 于 SIMLS/SILS 的 局 部 函数 和 相 
应 的 新 的 广义 有 限 元 插值 对 正弦 函数 的 拟 合 结 果 示 于 图 9。 从 图 中 看 到 ， 虽然 局 部 函数 本 质 
上 是 逼近 的 ， 但 导出 的 广义 有 限 元 却 是 全 局 插值 的 。 


小 结 : 移动 最 小 二 乘 本 1 
来 只 是 一 种 拟 合 有 逼近 ， 但 由 0.5 
此 导出 的 新 广义 有 限 元 却 是 
全 局 插值 的 ， 严 格 满足 克 罗 os 


e—e—e 
40 02 04\06 08 


0.5 
ERMEER Fa. Bah A 
最 小 二 和 具有 性 能 强健 (无 ee i 
RED ase CORN EDREENUEGAEO 7 (ES CFEM EM 
小 (如 SILS) Aes. 图 9， 使 用 SIMLS 局 部 允 近 的 广义 有 限 元 (GFEM) 


r=3h， 被 逼近 的 函数 为 正弦 函数 


关于 该 方法 更 详细 的 内 容 可 参考 文献 [46]。 本 文 主要 目的 是 以 该 方法 为 例 ， 抛 砖 引 玉 ， 
讨论 过 去 被 认为 代价 昂贵 的 一 类 算法 在 新 兴 众 核 处 理 器 上 可 能 出 现 的 新 的 特征 。 


5.2 数值 试验 
5.2.1 收敛 性 测试 
考虑 下 列 的 一 维 测试 问题 
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精确 解 取 为 
误差 的 最 小 二 乘 范 数 和 和 外 


收敛 测试 中 使 用 了 
LE 


LOG (lull) 


1 2 
LOG( DoFs ) 


(a) 最 小 二 


lu| - (f. (u^ (x)- «(x)faaJ. e 


五 个 一 致 加 密 的 规则 网 格 ， 分 别 为 20、40、80、 


元 网 格 。 单 元 片 大 小 取 r=ch, c=0... “r= 0h” 对 应 于 标准 


乘 范 数 
图 10. 拉 格 明日 型 局 部 通 近 


信 ， 


息 技术 快报 
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—-Au=f x 
u(0)=0 
Vu(1)=0 


u(x) =x(x- 


BIL BLY EXA F 


一 + 一 r=0h (有 限 元 ),2.0 
Or r=1h ( 拉 格 朗 日 ) ,2.4 
x r=2h ( 拉 格 朗 日 ) ,3.4 
m r=3h ( 拉 格 朗 日 ) ,4.4 
一 @ 一 r=4h ( 拉 格 朗 日 ) ,5.3 
o r=5h ( 拉 格 朗 日 ) ,6.3 
—e- r=10h ( 拉 格 朗 日 ) 

-- 参考 线 , 11.0 


* DoFs ， 自 由 度 


3 


= (fa(Vu" (x)- 


» 10. 


e (0,1) 


2 


1) e € 


Vu( faa] 
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(14) 


(15) 


(16) 


160 以 及 320 个 线性 


有 限 元 情形 。 


+ r=0h (有 限 元 ),1.0 
^ r=th ( 拉 格 朗 日 ) ,1.4 


r=2h ( 拉 格 朗 日 ) ,2.4 


1 r=3h ( 拉 格 朗 日 ) ,3.3 
-外 -r=4h ( 拉 格 朗 日 ) ,4.3 
-O- r=5h ( 拉 格 朗 日 ) ,5.3 
> r-10h ( 拉 格 朗 日 ) 

--- 参考 线 , 10.0 


* DoFs ， 自 由 度 


: —— 1=0h (有 限 元 ),2.0 

-2 ~ r=1h (simis) ,2.29 一 一 r=0h (有 限 元 ),1.0 
E a enemy 4n T -man mo) 214 
5.6 —- r-4h (simis) ,5.16 = 38t" r-3h (simis) ,3.19 
= -e~ r-5h (simis) ,6.05 iS. -Cr r-4h (simls) ,4.07 
e -8 * DoFs ， 自 由 度 9 4^ r-5h (simls) ,4.86 
40 * DoFs, ARE 

12 2 

LOG( DoFs ) LOG( DoF s ) 
(a) 最 小 二 乘 范 数 (b) 能 量 范 数 
图 11. SIMLS 型 局 部 逼近 《图 例 中 数字 为 收敛 率 ) 

单元 上 高 斯 积分 点 数 取 为 ct2。 收 和 敛 性 测试 结果 见 图 10 和 图 11， 分 别 对 应 于 两 种 不 同 
的 局 部 盘 近 构造 方式 。 

10 和 图 11 中 的 结果 显示 ， 单 元 片 大 小 每 增 大 1h， 新 广义 有 限 元 的 收敛 阶 通 常 都 会 
相应 提高 一 阶 ( 收 敛 速 率 见 图 例 )。 这 一 收敛 属性 类 似 于 现 有 广义 有 限 元 ， 不 同 点 在 于 : 新 
广义 有 限 元 在 网 格 尺寸 和 自由 度数 不 改变 的 情况 下 , 便 可 以 获得 不 同 阶 的 精度 和 收敛 性 。 这 
种 特征 过 去 多 见于 有 限 差分 方法 。 在 传统 的 有 限 元 或 现 有 广义 有 限 元 中 ， 要 提高 收敛 阶 ， 必 
须 在 相同 网 格 尺 寸 下 增加 自由 度 ， 或 在 单元 上 增加 新 的 节点 (传统 有 限 元 )， 或 增加 节点 上 
的 广义 自由 度 〈 现 有 广义 有 限 元 )。 


5.2.2 新 方法 在 新 兴 众 核 处 理 器 GPU 上 的 性 能 测试 
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超级 计算 面临 的 挑战 及 其 对 未 来 数值 算法 设计 的 可 能 影响 


由 于 总 体 方 程 求解 与 局 部 通 近 方式 无 关 , 下 面 我 们 主要 对 新 方法 形成 的 总 体 方 程 的 求解 
进行 测试 。 对 于 单元 的 计算 与 组 装 ， 由 于 其 与 局 部 通 近 方式 有 关 ， 和 暂 不 考虑 〈 但 不 影响 结果 
的 一 般 性 )。 

首先 考察 新 方法 与 现 有 广义 有 限 元 方法 在 总 体 刚 度 矩 阵 形式 上 的 差别 。 从 图 12 中 对 新 
方法 与 现 有 方法 总 体 刚 度 和 矩阵 特征 的 比较 ， 我 们 观察 到 : 


- ”在 新 方法 中 ， 总 体 和 矩阵 的 大 小 不 随 局 部 盟 近 阶 次 的 改变 而 改变 ， 而 在 现 有 方法 中 ， 
总 体 和 矩阵 的 维 数 会 随 着 局 部 交 近 阶 次 的 提高 而 增 大 。 

- ”在 新 方法 中 ， 总 体 矩 阵 的 稀 琉 性 〈 非 零 元 素 占 比 ) 随 着 局 部 远近 阶 次 的 提高 而 降 
低 ， 而 在 现 有 方法 中 ， 改 变局 部 冯 近 的 阶 次 无 法 改变 总 体 和 矩阵 的 稀 玉 C 性 。 


WY 


(a)FEM(25.6%) (b)r-h(53.796) (c)r-2h(75.296) (d)r=3h(90.0%) 
(e)FEM(25.6%)  (f)p-1(12.896) (g)p=2(12.5%) (h)p=3(12.1%) 


图 12. 一 维 问题 10 E RT INE BER s Hc SERERE BE 
(a)-(d) AMA, (e)-(h) AMA AIK. 括 绝 中 的 数字 为 非 替 元 素 占 比 , RORY HR 
性 . DUA A IK AY SAREE A SSE YOR AR m va JU 78 RARE JR XA, 而 在 新 方 
法 中 ,， 随 着 局 部 逼近 阶 次 的 提高 , 总 体 矩 阵 变 得 稠密 . 


在 新 方法 中 ， 总 体 矩 阵 的 稀 玻 性 是 变化 的 或 可 调节 的 。 这 也 是 目前 无 网 格 /粒子 类 方法 
的 一 个 共同 特征 。 因 此 ， 下 面 的 结论 对 于 无 网 格 /粒子 类 方法 具有 一 定 的 普遍 意义 。 


测试 问题 为 二 维 泊 松 问题 。 测试 平台 为 GPU 众 核 和 CPU 单 核 (只 用 一 个 核 进行 计算 以 
反映 PC 时 代 的 计算 特征 )。 性 能 指标 为 单位 时 间 内 浮 点 运算 次 数 Gflop/s。 测 试 对 象 为 总 体 
方程 组 的 求解 环节 〈 使 得 测试 更 具有 一 般 性 )。 求 解 方法 采用 共 斩 梯 度 向 量 法 。 求 解 器 采用 
基于 CUDA 的 开源 线性 代数 库 CUSP。 稀 疏 矩 阵 采 用 CSR 存储 格式 。 测 试 结果 示 于 图 13. 


首先 ， 令 r=0h 将 新 方法 退化 为 标准 有 限 元 方法 ， 对 标准 有 限 元 方法 在 CPU 单 核 (PC 
时 代 ) 与 GPU 众 核 上 的 性 能 进行 比较 , 结果 见 图 13 中 圆圈 中 的 两 个 数据 点 。 从 图 中 我 们 观 
察 到 , 虽然 CPU 单 核 与 GPU 众 核 的 峰值 性 能 相差 很 大 (前 者 为 10Gflop/s, 后 者 为 78Gflop/s )， 
但 标准 有 限 元 在 两 种 计算 平台 上 的 性 能 却 相差 无 儿 。 这 一 测试 结果 表明 : 如 果 在 众 核 上 使 用 
标准 有 限 元 算法 ， 众 核 强大 的 计算 能 力 将 可 能 得 不 到 充分 发 挥 和 利用 《算法 计算 密集 度 低 ， 
性 能 受 访 存 限制 )。 


AS 


ce 
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PR, RITEAR One PATE Ce I A), HII BR te PE AE 
CPU 单 核 和 GPU 众 核 上 重新 进行 性 能 比较 (图 中 实 线 )。 我 们 发 现 ， 在 GPU 众 核 上 ， 当 提 
高 方法 的 阶 次 ， 即 增 大 单元 片 大 小 了 《或 降低 总 体 窍 阵 的 稀 玻 性 )， 新 方法 的 性 能 得 到 非常 
明显 的 提升 “硬件 性 能 得 到 了 利用 和 发 挥 )， 与 方法 阶 次 几乎 呈 线 性 增长 。 这 一 测试 结果 表 
BH: 新 方法 的 高 阶 格式 能 充分 发 挥 众 核 的 计算 能 力 。 然 而 ， 图 12 显示 新 方法 导致 较为 稠密 
的 总 体 和 矩阵 ， 如 果 从 PC 时 代 的 传统 观点 看 ， 这 类 方法 应 该 是 较为 低 效 的 算法 。 上 述 测试 却 
表明 ， 在 众 核 上 由 于 硬件 性 能 得 到 了 充分 利用 ， 这 类 方法 反而 变 得 更 加 高 效 。 


图 13 中 虚线 表示 同样 的 高 价 方法 在 _ 14 


—m— 众 核 GPU 


CPU 单 核 上 的 测试 结果 。 我 们 看 到 , 如 果 在 B12 o- 单 核 CPU 

CPU 单 核 上 (或 PC 时 代 )， 高 阶 格式 并 不 O10 NEM = 
4 x. r= 

能 提高 硬件 的 计算 性 能 利用 率 , 反而 会 导致 m dá e 


By MESS Mr ee aa Ai pos] 
du E ert 
代 , FE AT BO eS BOAT ETT IE OR r=ch, c>0 


认为 是 低 效 或 计算 费力 的 事实 是 相符 的 (在 p) des ee ee a 


PC/N Bi e, FECCE IY EEO EEE Nr (2.3. A S S 
单 核 时 代 ， 硬 件 的 计算 能 力 是 主要 瓶 oy a a Re 


由 于 新 广义 有 限 元 方法 ， 无 网 格 方法 ， 
以 及 粒子 方法 等 一 类 新 兴 计 算 方法 (外 加 有 BB. 新 方法 在 众 
限 差分 方法 ) 具有 在 不 改变 网 格 的 前 提 下 便 
可 获得 不 同 精度 的 共性 特点 〈 其 总 体 和 矩阵 均 具 有 图 12(9)-(d) 的 变化 特征 )。 我 们 相信 ， 前 述 
实验 观察 到 的 新 广义 有 限 元 方法 在 众 核 处 理 器 上 的 计算 特性 极 有 可 能 具有 一定 的 普 适 性 。 


核 平台 上 的 性 能 测试 


依据 PC 时 代 的 经 验 ,， 我 们 总 是 假设 浮 点 运算 能 力 是 宝贵 的 ， 所 以 习惯 于 以 频繁 访 存 来 
避免 重复 计算 , 算法 设计 和 优化 的 主要 原则 也 是 减少 浮 点 运算 次 数 。 但是， 从 目前 的 技术 发 
展 趋势 看 , 未 来 科学 与 工程 应 用 所 使 用 的 计算 平台 的 硬件 计算 能 力 与 访 存 能 力 的 比率 必 将 越 
来 越 高 。 因 此 从 硬件 角度 看 ,设计 算法 时 单位 访 存 所 进行 的 浮 点 运算 越 多 越 好 。 这 与 过 去 科 
学 与 工程 应 用 算法 设计 的 思考 方向 和 编程 习惯 正好 相反 。 对 于 领域 应 用 专家 而 言 ， 有 必要 结 
合 人 硬件 的 特征 和 变化 趋势 ， 对 关键 数值 算法 进行 重新 审视 、 选 优 ， 甚 至 重新 设计 ; 在 传统 的 
精度 、 收 敛 性 等 指标 之 外 ， 重 视 计算 密集 度 、 数 据 移 动 复杂 性 等 算法 特征 。 


本 文 抛砖引玉 ， 面 向 新 兴 计 算 机 体系 结构 ， 从 充分 释放 众 核 处 理 器 “ 见 余 ” 计 算 能 力 的 角 
度 出 发 ， 探 索 了 针对 下 一 代 处 理 器 结构 发 展 高 效 数值 算法 的 一 个 可 能 的 技术 思路 。 
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